GRRData<-function(x,k,n,r){
    #Create by Phoebe.Chang<phoebe.chang@intelligentgroup.cn> 2016-06-14
    #处理空数据--By Phoebe 2016-07-13
    x[is.nan(x)]<-NA
    # 按维度组成数据框
    testData<-data.frame(testResult=x[1:(k*n*r)],
                            Appraisers=rep(gl(k,r,k*r),n),
                            Parts=gl(n,k*r,k*n*r))
  return(testData)
}

DataValidity<-function(x,k,n,r){
    #Create by Phoebe.Chang<phoebe.chang@intelligentgroup.cn> 2016-06-29
    # 按维度组成数据框
    testData<-GRRData(x,k,n,r)
    # 按人维度判定是否有一个维度全为NaN
    aV<-aggregate(testData$testResult,list(testData$Appraisers),mean,na.rm=T)$x
    # 按产品维度判定是否有一个维度全为NaN
    pV<-aggregate(testData$testResult,list(testData$Parts),mean,na.rm=T)$x
    if(sum(is.nan(aV))==0 & sum(is.nan(pV))==0 & (k>1 | n>1)){
        #两个因素维度不能同时为1某一个维度不能全是空值
        Validity<-TRUE
    }else{
        Validity<-FALSE
    }
    return(Validity)
}

GRRAnova<-function (x,k,n,r){
    #Create by Phoebe.Chang<phoebe.chang@intelligentgroup.cn> 2016-05-16
    # 设置数据精度
    options(digits=10)
    # 按维度组成数据框
    testData<-GRRData(x,k,n,r)
    #判断数据有效性--2016-06-29
    validity<-DataValidity(x,k,n,r)
    if(validity){
    #根据不同情况确定分析模型
        if(k<2){
            #自动测试，关于产品的单因分析
            fit<-aov(testResult~Parts,data=testData)
            #统一输出格式
            anovaData<-summary(fit)[[1]]
            #只测试1次的情况,统一输出格式
            if(r<2){
                anovaData[,"F value"]<-rep(NA,1)
                anovaData[,"Pr(>F)"]<-rep(NA,1)
                anovaData["Residuals",]<-rep(NA,5)
            }
            add<-rep(NA,length(anovaData))
            anovaData<-rbind("Appraisers"=add,anovaData[1,],"Appraisers:Parts"=add,anovaData[2,])
        }else{
            if(n<2){
                #无产品的GRR测试，关于人的单因分析
                fit<-aov(testResult~Appraisers,data=testData)
                #统一输出格式
                anovaData<-summary(fit)[[1]]
                add<-rep(NA,length(anovaData))
                anovaData<-rbind(anovaData[1,],"Parts"=add,"Appraisers:Parts"=add,anovaData[2,])
            }else{
                # 双因交叉方差分析，分析模型Appraisers + Parts + Parts:Parts
                fit<-aov(testResult~Appraisers+Parts+Appraisers:Parts,data=testData)
                anovaData<-summary(fit)[[1]]
                #只测试1次的情况,统一输出格式
                if(r<2){
                    anovaData[,"F value"]<-rep(NA,3)
                    anovaData[,"Pr(>F)"]<-rep(NA,3)
                    anovaData["Residuals",]<-rep(NA,5)
                }
            }
        }
    }else{
        #当数据不合法时组合成空数据框
        Df<-c(k-1,n-1,(k-1)*(n-1),n*k*(r-1))
        x1<-rep(NA,4)
        anovaData<-data.frame(Df)
        row.names(anovaData)<-c("Appraisers","Parts","Appraisers:Parts",
                                "Residuals")
        anovaData[,"Sum Sq"]<-rep(NA,4)
        anovaData[,"Mean Sq"]<-rep(NA,4)
        anovaData[,"F value"]<-rep(NA,4)
        anovaData[,"Pr(>F)"]<-rep(NA,4)
    }
    #跟着Minitab一起错让FA=MSSA/MSSAP、FB=MSSB/MSSAP--by Phoebe
    anovaData$`F value`[1]<-anovaData$`Mean Sq`[1]/anovaData$`Mean Sq`[3]
    anovaData$`F value`[2]<-anovaData$`Mean Sq`[2]/anovaData$`Mean Sq`[3]
    anovaData$`Pr(>F)`[1]<-1-pf(anovaData$`F value`[1],anovaData$Df[1],anovaData$Df[3])
    anovaData$`Pr(>F)`[2]<-1-pf(anovaData$`F value`[2],anovaData$Df[2],anovaData$Df[3])
    #判断交叉因素影响是否显著
    p<-anovaData$"Pr(>F)"[3]
    if(p>0.05 & !is.na(p)){
    #交叉因素影响不显著，不考虑交叉因素影响
    anovaData$"Mean Sq"[4]<-(anovaData[3,2]+anovaData[4,2])/(anovaData[3,1]+anovaData[4,1])
    }
    # 获取并加工方差分析结果
    len<-length(anovaData)-2
    temp<-c(sum(anovaData[,1][!is.na(anovaData[,1])]),
            sum(anovaData[,2][!is.na(anovaData[,2])]),
            rep(NA,len))
    #return(temp)
    anovaData["Total",]<-temp
    # 返回方差分析计算结果
    return(anovaData)
}

GRRSource<-function (x,k,n,r,sig=6,tole){
    #Create by Phoebe.Chang<phoebe.chang@intelligentgroup.cn> 2016-05-16
    #Last Update by Phoebe.Chang<phoebe.chang@intelligentgroup.cn> 2016-06-22
    #获取方差分析结果
    an<-GRRAnova(x,k,n,r)
    # 获取Appraiser*Part的检验P值
    PAP<-an$"Pr(>F)"[3]
    #获取均方差
    MS<-an$"Mean Sq"
    #计算出方差序列
    if(PAP>0.05 | is.na(PAP)){
        #当P值>0.05,或者P值为空接受A*P交叉因素对结果影响不大的假设不考虑A*P交叉因素——By Phoebe 2016-06-22
        va<-c(MS[4],NA,(MS[1]-MS[4])/(n*r),NA,NA,(MS[2]-MS[4])/(k*r),NA)
    }else{
        va<-c(MS[4],NA,(MS[1]-MS[3])/(n*r),(MS[3]-MS[4])/r,NA,(MS[2]-MS[3])/(k*r),NA)
    }
    #处理负数--Fix by Phoebe<2016-05-27>
    va[va<0]<-NA

    va[2]<-sum(va[3:4][!is.na(va[3:4])])
    va[5]<-sum(va[1:2][!is.na(va[1:2])])
    va[7]<-sum(va[5:6][!is.na(va[5:6])])
    #计算出Sigma
    Sigma<-sqrt(va)
    #计算nSigma
    nSigma<-sig*Sigma
    #计算总变差比率Contribution未加%
    Contribution<-va/va[7]*100
    #计算总方差比率Variance未加%
    Variance<-Sigma/Sigma[7]*100
    #计算贡献率Tolerance未加%
    Tolerance<-nSigma/tole*100
    #Deal with the negative ——fix by phoebe<2016-11-04>
    Tolerance[Tolerance<0|Tolerance==Inf]<-NA
    GRRData<-data.frame(va,Sigma,nSigma,Contribution,Variance,Tolerance)
    row.names(GRRData)<-c("Repeatability","Reproducibility","Appraisers",
                            "Appraisers x Part","Gage R&R","Parts","Total")
    return(GRRData)
}

GRRXbar<-function(x,k,n,r){
    #Create by Phoebe.Chang<phoebe.chang@intelligentgroup.cn> 2016-06-13
    #Last Update by Phoebe.Chang<phoebe.chang@intelligentgroup.cn> 2016-06-14
    # 按维度组成数据框
    testData<-GRRData(x,k,n,r)
    # 分组计算均值
    Xbar<-aggregate(testData$testResult,list(testData$Parts,testData$Appraisers),mean,na.rm=T)
    # 获取结果数组
    Xbar<-Xbar$x
    # 保存参数表--Phoebe 2016/06/14
    A2<-c(1.880,1.023,0.729,0.577,0.483,0.419,0.373,0.337,0.308,0.285,0.266,0.249,
        0.235,0.223)
    # 计算上中下控线--Phoebe 2016/06/14
    if(r<2|r>15){
        UCL<-NA
        LCL<-NA
    }else{
        Rdata<-GRRRange(x,k,n,r)[k*n+2]
        UCL<-mean(Xbar)+Rdata*A2[r-1]
        LCL<-mean(Xbar)-Rdata*A2[r-1]
    }
    # 组合计算结果--Phoebe 2016/06/14
    Xbar<-c(Xbar,UCL,mean(Xbar),LCL)
    return(Xbar)
}

GRRRange<-function(x,k,n,r){
    #Create by Phoebe.Chang<phoebe.chang@intelligentgroup.cn> 2016-06-13
    #Last Update by Phoebe.Chang<phoebe.chang@intelligentgroup.cn> 2016-06-14
    # 按维度组成数据框
    testData<-GRRData(x,k,n,r)
    # 定义极差函数
    Range<-function(x, na.rm = FALSE){
        #处理空数据--By Phoebe 2016-07-13
        if(na.rm){
            x <- x[!is.na(x)]
        }
        diff(range(x))
    }
    # 分组计算极差
    Rdata<-aggregate(testData$testResult,list(testData$Parts,testData$Appraisers),Range,na.rm=T)
    # 获取结果数组
    Rdata<-Rdata$x
    # 保存参数表--Phoebe 2016/06/14
    D4<-c(3.267,2.575,2.282,2.115,2.004,1.924,1.864,1.816,1.777,1.744,1.716,1.692,
            1.671,1.652)
    # 计算上中下控线--Phoebe 2016/06/14
    if(r<2|r>15){
        UCL<-NA
    }else{
        UCL<-mean(Rdata)*D4[r-1]
    }
    # 组合计算结果--Phoebe 2016/06/14
    Rdata<-c(Rdata,UCL,mean(Rdata),0)
    return(Rdata)
}

GRRAppraiser<-function(x,k,n,r){
    #Create by Phoebe.Chang<phoebe.chang@intelligentgroup.cn> 2016-06-15
    # 按维度组成数据框
    testData<-GRRData(x,k,n,r)
    # 按人排序
    appraiser<-testData[order(testData$Appraisers,decreasing = F),]
    #形成数据矩阵
    appraiserData<- matrix(appraiser$testResult, nrow=k,ncol= n*r,byrow=TRUE)
    #计算均值矩阵
    ave<-matrix(apply(appraiserData,1,mean,na.rm=T))
    #合并矩阵
    appraiserResult<-cbind(appraiserData,ave)
    return(appraiserResult)
}

GRRPart<-function(x,k,n,r){
    #Create by Phoebe.Chang<phoebe.chang@intelligentgroup.cn> 2016-06-15
    # 按维度组成数据框
    testData<-GRRData(x,k,n,r)
    # 按产品排序
    part<-testData[order(testData$Parts,decreasing = F),]
    #形成数据矩阵
    partData<- matrix(part$testResult, nrow=n,ncol= k*r,byrow=TRUE)
    #计算均值矩阵
    ave<-matrix(apply(partData,1,mean,na.rm=T))
    #合并矩阵
    partResult<-cbind(partData,ave)
    return(partResult)
}

GRRPAPlot<-function(x,k,n,r){
    #Create by Phoebe.Chang<phoebe.chang@intelligentgroup.cn> 2016-06-22
    # 按维度组成数据框
    testData<-GRRData(x,k,n,r)
    #分组计算均值
    Xbar<-aggregate(testData$testResult,list(testData$Parts,testData$Appraisers),mean,na.rm=T)
    #按维度生成矩阵
    PAPlot<-matrix(Xbar$x,nrow=k,ncol=n,byrow = TRUE)
    return(PAPlot)
}

GRRComponent<-function (x,k,n,r,sig=6,tole){
    #Create by Phoebe.Chang<phoebe.chang@intelligentgroup.cn> 2016-06-22
    # 获取Source表的结果
    sourceData<-GRRSource(x,k,n,r,sig,tole)
    # 去掉不需要的数据
    sourceData<-sourceData[c(-3,-4,-7),]
    sourceData<-sourceData[,c(-1:-3)]
    # 转化为矩阵
    componentResult<-t(as.matrix(sourceData))
    rownames(componentResult)<-paste('%',rownames(componentResult))
    return(componentResult)
}

